Linearized Auxiliary fields Monte Carlo: efficient sampling of the fermion sign 
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We introduce a method that combines the power of both the lattice Green function Monte Carlo 
(LGFMC) with the auxiliary field techniques (AFQMC), and allows us to compute exact ground 
state properties of the Hubbard model for U < 4t on finite clusters. Thanks to LGFMC one 
obtains unbiased zero temperature results, not affected by the so called Trotter approximation of 
the imaginary time propagator e~^^. On the other hand the AFQMC formalism yields a remarkably 
fast convergence in r before the fermion sign problem becomes prohibitive. As a first application 
we report ground state energies in the Hubbard model at U/t = 4 with up to one hundred sites. 
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After several years of scientific effort, based on ad- 
vanced analitycal and numerical methods, only very few 
properties of the 2D Hubbard model have been settled. 
The 2D Hubbard model is defined in a square lattice con- 
taining a finite number L (N) of sites (electrons): 

H=-t 4,.c,,a + Uj2nJni (1) 

<i,j>,(7 i 

with standard notations. In the thermodinamic limit, 
namely for L — >■ cx) at given density p = N/L funda- 
mental issues such as the existence of a ferromagnetic 
phase at large U/t ratio and/or the stability of an homo- 
geneous ground state with possible d-wave supoercon- 
ducting properties are still highly debated, as several ap- 
proximate numerical techniques lead to controversial and 
often conflicting results. This situation is particularly 
embarazzing, since recent progress in the realization of 
fermionic optical lattices could lead soon to the experi- 
mental realization of the fermionic Hubbard model, ap- 
parently much before we will reach a consensus among 
the different theoretical and numerical techniques. 

Method: In LGFMCi, the main property used to com- 
pute the ground state tjjo of a many body Hamiltonian H, 
is the iterative application of a linearized Green's func- 
tion: 

G = KI -H (2) 

to an initial wave function -0 by a stochastic method, 
namely ■f/^o — lim G'"|V')- Here / is the identity matrix 

and A is a suitably large constant. This is possible be- 
cause the application of G to a given configuration \x), 
where electrons have definite positions and spins, can be 
expressed as a sum of a finite number of independent 
configurations. Namely the number of non zero matrix 
elements Gx' ,x — {x'\G\x) for given x is affordable (cx i), 
though the Hilbert space is exponentially large with L. 
For large n it is possible to sample the ground state wave 
function 'i/'ol^;) = (s^lV'o) and its correlation functions. 
This method can be improved by the so called impor- 
tance sampling, yielding a much more efficient algorithm: 



the matrix elements of G are scaled by means of the so 
called guiding function ipg{x): 

Gx',x ^ ■ipgix')Gx',x/lpg{x) (3) 

This method allows the calculation of exact ground state 
properties without any approximation other than (i) fi- 
nite L and (ii) statistical errors. The latter may be partic- 
ularly large with fermionic systems, due to the unfamous 
"sign problem", a limitation that is usually much severe, 
if not prohibitive, for this type of approach. We remark 
also that, in LGFMC, one can work with infinite and 
sample the many body propagator e~^^ — lim (x)" 

with a similar computational effort, despite n — rA — > oo 
in this limit. Moreover, since the projection time r is pro- 
portional to the length of the Monte Carlo simulation, 
converged ground state properties are easily obtained af- 
ter the equilibration time, and the limit r — cx) is basi- 
cally achieved without any particular effort when there 
is no sign problem. 

Another important stochastic method that is quite 
popular for the Hubbard model is based on AFQMC— 
and the algebra of one body propagators. The basic prop- 
erty used in this approach is that a one body operator 
U, when applied to a Slater determinant \D), generates 
again a Slater determinant \D') — U\D), that is easy to 
evaluate. Obviously, by the above definition, the prod- 
uct of several one-body propagators remains a one-body 
propagator, and the computation is always feasible for 
a finite number of them. Indeed, within AFQMC, the 
many body propagator e~^'^ is conveniently written as 
a superposition of time dependent one body propagators 
U[a]{T,0) = Ua{tL^) ■ ■ ■ Ucr{ti) after the introduction of 
discrete Hubbard- Stratonovich time-dependent auxiliary 
fields (Ji{tj) = ±1 defined in all the L lattice sites and 
a finite number of discrete imaginary times tj — jr / L-r 
j ^ I,-- - ,Lr. Thus e-"^\i;) = Y.UaiT,0)\tp) for large 

imaginary time gives formally the exact ground state 
wave function and can be computed by Monte Carlo sam- 
pling of the auxiliary fields cri{tj) = ±1— At variance of 
the LGFMC technique, it is not simple to avoid the er- 
ror due to the discretization in time of the propagator- 
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usually called Trotter error- and all the results require a 
careful and often boring extrapolation both in r — >■ c» 
and Lr ^ CO. By contrast, within AFQMC, the sign 
problem is much less severe for two main reasons, (i) for 
U = the method is exact and one does not need any 
random sampling, (ii) at half-filling or with negative U 
the method is not affected by the " sign problem" . 

In the following we propose to combine the power of the 
two methods, by taking the best of the two approaches 
in what we name Linearized Auxiliary fields Monte Carlo 
(LAQMC). From AFQMC we take the advantage of a 
much less severe sign problem and from LGFMC the ex- 
act imaginary time projection will be available in a rig- 
orous and simple way. The latter achievement has been 
made recently possible also within the so called "dia- 
grammatic" Monte Carlo§'-', but only within the reason- 
able assumption that the diagrammatic perturbation se- 
rie converges. 

In order to define this new approach we use that the 
lattice Green function © of the Hubbard model can be 
written exactly as a finite sum of one body propagators 
Ui, with an approach that is similar to the conventional 
auxiliary field technique, where instead of splitting up the 
many propagator exp(— tTJ) we focus on its linearized 
expression given by G: 



G = hi - H = ^a,U, 



(4) 



where the coefficients ai > and the Ui will be de- 
fined later on and p = 4L for the Hubbard model. The 
above identity can be generally fuUfiUed for any reason- 
able many-body Hamiltonian by using a number p of 
one-body propagators that scales at most as L^. For 
the Hubbard model we obtain after simple algebra: 

p 

A = Tuie^ + e-^ ^2)N + ^a^ 

1=1 



U, 



+5/2L±L/2 



where, in order to satisfy Eq.Q: 

a,(e-<'-l) = 
2a,(e^-l)(e-^-l) = 



-U i>2L 



(5) 
(6) 



and ef.i ~ — 2t(cos(fc^) -|- cos(fc^)) is the U — band and 
fc* label all the 2L independent k vectors of the spin- 
up and spin-down electrons {ai =t for i — 1 , • • • i and 
(Tj =1 for i = L + 1, • • • 2L). Here for simplicity we set 
ai = Tt {ai — Tu), independent of i for i < 2L (i > 
2L). Notice that the main difference compared with the 
discrete Hubbard- Stratonovich transformation^ is that, 
in order to decompose the propagator G, we introduce 
not only one body propagators for the interaction term 
U, but we use also further one body operators Ui for 
i < 2L, to recast the kinetic term as a simple sum of 



one-body propagators. In some sense this is equivalent to 
double the dimension of the auxiliary fields [a] , extension 
that does not lead to a particular loss of efficiency of the 
algorithm and, on the other hand, allows us to remove 
the bias due to the time discretization in a simple and 
rigorous way. The choice of Tt and Tu , and in principle 
all the coefficients a^, are completely arbitrary in this 
approach and can be tuned for optimizing efficiency, by 
a substantial alleviation of the sign problem. On the 
other hand, it is simple to realize that for Tt = 1/At 
and A = VUAt, one recover the same Master equation^ 
of the standard AFQMC in the limit Ai — > 0, where 
At is the time discretization adopted with the Trotter 
approximation. Thus at least in this limit the proposed 
method has no sign problem in all the cases when the 
standard AFQMC has no sign problem. We will refer 
as the "time continuous limit" (TCL) for this particular 
choice of the coupling Tt and A. 

Once this decomposition is implemented we immedi- 
ately recover the main property of GFMC, provided the 
configuration \x) is replaced by a Slater determinant \D), 
defined by the orbital matrix Dij, that in the particular 
case of the Hubbard model reads: 



N 



i=l j i=JVT + l ] 

where D is real and is the number of spin- up particles. 
In fact, the application of G, written in the form to 
a single Slater determinant generates a finite number of 
Slater determinants of the same form. The importance 
sampling can be analogously defined by means of a guid- 
ing function : D ^ ^g{D) such that ipg{D) can be 
easily computed. From this point of view the method is 
similar to Constrained Path AFQMCS, where the guid- 
ing function is defined in terms of a simple mean field 
Slater determinant IV'mf), namely ipg{D) = (^Jjmf\D)- 

In order to work with a rigorous statistical method 
with finite variance one can employ standard smearing 
procedures of the guiding function based on the reweight- 
ing method, that allow us to work with a non zero guid- 
ing function for all D generated by the Markov process^. 
Following the argument discussed in RefjlO), a very effi- 
cient smearing procedure is defined here by means of the 
Green's function, a 2L x 2L matrix: 



{'>pMF\ci^Cj^^'\D) 

{iPmf\D) 



(8) 



If the determinant |D| — {iPmf\D vanishes, some of the 
the Green's function elements should diverge in the same 
way, so that we define a reweighting factor i? < 1, satis- 
fying R ~ \D\ for \D\ -> 0: 



R{D) = l/Jl + e^TriggT) 



(9) 



where the trace Tr{gg^) = J2i,j,a,cr' Ifia.ja'P, and e is 
an appropriate small number, chosen in a way that, on 
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FIG. 1: (Colors online) Comparison of exact ground state 
energy and the one obtained by LAQMC projection for the 
Hubbard with two holes (A'^ = 16) in the 18 square lattice with 
Ft = 2.5(7. Notice the rapid convergence with projection time 
even for large U/t. 



average, R is around 0.9 0.95. Then it is easy to de- 
fine a guiding function ipg{D) — ^ that remains finite 
for \D\ — T' 0. On the other hand for e small enough -0^ 
remains sufficiently close to the mean field determinant 
|D|, by allowing a very good importance sampling. 

Let us consider the basic step of the stochastic imple- 
mentation of the power method A — _ff . Once the nor- 
malization bo = X^iLi l^;7i|-D).|£i) I is given, the walker 
weight and the determinant \D) are updated by means 
of the following simple Markov chain: 

\D) ^ Ui\D) with probability \Gui\D),\D)\lbD 
w wbDSgnGui\D),\D) (10) 

On the other hand the expectation value of the energy 
as any other " mixed average" quantity, can be computed 
by means of the ratio of two random variable averages 
'^"'^krI^^ , where eL{D) is the local energy: eL{D) = 

Several walkers with weights Wi,i = 1,---M evolve 
with the above Markov chain and undergo a " branching" 
processii to optimize efficiency in the sampling. In order 
to eliminate completely the finite population bias it is 
necessary only to bookkeep a "correcting factor" w — 

Constrained path as a standard LGFMC: In this 
formalism it is very simple to employ the constrained 
path approximation^. This is a very stable algorithm as 
the walker is constrained to avoid regions with extremely 
small determinant \D\. A simple way to implement the 
CPQMC approximation in the continuous limit is ob- 
tained by the following standard recipe: whenever a sign 
change occur in Eq. ljlOl) . the walker weight w is simply 
annihilated within this approximation. This implemen- 
tation is particularly important for applying the standard 
release nodc^, explained below. 

Release node technique: When there is sign prob- 
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FIG. 2; (Colors online) Energy as a function of the power 
method iterations using different mean field wave functions 
'4^MF for defining the guiding function in LAQMC. They differ 
for the value of the chemical potential and the BCS d- 
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wave parameter A^Qg (see RefllSl). (a) /io ~ — 0.198t (b) 
Agcf = 0.002t 



lem, the Markov chain ([TU)) is unstable because after a 
while half of the walkers will have negative sign and will 
cancel almost exactly the contribution of the ones with 
positive sign. In order to stabilize the process, we em- 
ploy the standard release node technique introduced long 
time ago^, and adapted to the present case, to take into 
account only a "discrete time" projection given by the 
power method. Therefore after the Markov chain equili- 
brates, we can have access with a single run and a very 
simple postprocessing, to all the history evolution of the 
energy as a function of the power method iterations start- 
ing from the very good estimate provided by the CPQMC 
state ipcPQMC, namely: 



En 



{lpMF\HC\ijCPQMc) 
{lpMF\G"-\lpCPQMc) 



(11) 



for all n < rir, where 71^ is the maximum release time 
allowed. 

Guiding function: In order to define the appropriate 
guiding function ipg at one electron per site filling N = L 
we use an antiferromagnetic mean field Slater determi- 
nant with the order parameter along the x-spin axisi^ii^. 



Moreover away from half-filling we use also a BCS 
like wave function with a small d-wave order parame- 

2 _ 2 

ter A^g(jg , when the ground state of the U — Q model 
is degenerate (open shell case). In the half-filled case 
iV = L in any bipartite lattice and for [/ = 0, there is no 
sign problem (for Ft large enough) and i^cPQMC = "00 ■ 
Therefore the method is exact already for n = 0. On 
the contrary the TCL will lead to biased result without 
smearing the guiding function (e — 0), due to lack of 
ergodicity in the Markov process. This is a well known 
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FIG. 3: (Colors online) Energy per hole for LAQMC release 
nodes (full symbols) and CPQMC (empty symbols). Lines 
are guides to the eye. 



problem in the standard CPQMC—, that can be defini- 
tively solved by means of the careful regularization in- 
troduced for e > 0. As a particular example we show in 
Fig.([T]) the comparison of the exact results^^, for the two 
hole case in the 18 site cluster. In this picture we reach 
convergence within statistical errors always with a small 
number of power iterations with no particular difficulty 
to sample the sign even for large U /t values. From this 
picture we remark that only for very large U jt ratios we 
observe that the convergence to the exact result is non 
monotonic as a function of n. In Fig. ^ the method 
is shown to work quite effectively as the CPQMC esti- 
mate remains very accurate also for a large cluster size, 
and convergence can be achieved much more quickly than 
the standard AFQMC^**. In this figure we see that even 
when the initial wave function is not optimal it is pos- 
sible to reach convergence with a quite good error bar. 
However it is also clear that a good initial guess allows 
a much more accurate energy estimate by using a much 
smaller rir- 



Finally we show the results^S, obtained for the energy 
per hole in Fig. ([3]), where ^ = 1 — N jL is the hole con- 
centration. If there is phase separation between an hole 
rich phase and the undoped insulator at (5 = the en- 
ergy per hole should have a minimum at a critical doping 
i5i~. Our reference energy at (5 = is exact and given by 
EjL = -0.85996(5)* for L ^ oo at U jt = 4. For smafi 
doping, though we obtained a non monotonic behavior 
similar to the ones shown in Fig. [1] for U jt > 12t, we 
were able to achieve convergence even in these particu- 
larly difficult cases using up to Ut — 10000 power method 
iterations. 

Considering therefore our energy results for the larger 
possible cluster size, we find a clear flat region in the 
energy per hole eh{S), as shown in Fig. ([3]), that may 
suggest an incipient phase separation^^. This behavior 
is in disagreement with the standard CPQMC, where 
a clear minimum was found at doping Sc — 5%^ for 
L ~ 100. In order to understand this result, we have 
performed CPQMC calculations with our technique, and 
found, as clearly displayed in the same picture, that the 
CPQMC energy per hole is quite sensitive to the quality 
of the guiding function in the small doping region, and a 
magnetic guiding function containing both an antiferro- 
magnetic and a d-wave order parameter greatly improves 
the CPQMC energy per hole, that becomes qualitatively 
similar to the LAQMC eh{S). Summirizing an infinite 
compressibility appears plausible only when the doping 
approaches zero, in contrast with the clearly bounded 
one obtained in the t — J model with a similar guiding 
functioi*i^. We believe that this difference is due to the 
particle-hole symmetry, that is not satisfied in the t — J 
model, and instead could imply eh{5) ~ 5^ as in ID. 

In conclusion we have presented a new method for the 
simulation of the Hubbard model with about 100 elec- 
trons, even when the sign problem prevents affordable 
calculations with standard techniques. Preliminary re- 
sults show a very large but finite compressibility at small 
doping in the square lattice Hubbard model at U /t — 4. 
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